function results = osl_get_recon_timecourses( source_recon_results, mask_fname, noreject,mni_coords )
   
% results = osl_get_recon_timecourses( source_recon_results )
%
% results = osl_get_recon_timecourses( source_recon_results, mask_fname )
%
% returns results.source_timecourses (nvox x ntrials x ntpts)
%
% MWW 2012

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% set up mask indices from mask/coords

results=[];
results.source_recon=source_recon_results.source_recon;
source_recon_name='source_recon';

if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis

    error('Incompatible with working in sensor space');

else

    if exist('mni_coords'),
        results.mask_indices_in_source_recon=1:size(mni_coords,1);
        results.mni_coord=mni_coords;
    else
        % setup std space brain
        if isfield(source_recon_results,'gridstep') && ~isempty(source_recon_results.gridstep)
        results.gridstep=source_recon_results.gridstep;

        S=[];
        S.lower_level_mask_fname=[results.source_recon.dirname '/' source_recon_name '_mask'];
        S.current_level_mask_fname=[results.source_recon.dirname '/' 'get_recon_timecourse_mask'];

        if ~exist('mask_fname','var') || isempty(mask_fname)
            mask_fname=S.lower_level_mask_fname;
        end;

        S.current_level.mask_fname=mask_fname;
        S.lower_level_mni_coord=source_recon_results.mni_coord;
        S.lower_level_gridstep=source_recon_results.gridstep;

        %S.current_level=first_level;

        [results.mask_indices_in_source_recon, results.mni_coord]=setup_mask_indices(S);

        results.mask_fname=S.current_level.mask_fname;

        clear S;

        elseif isfield(source_recon_results,'mni_coord') && ~isempty(source_recon_results.mni_coord)
          results.mni_coord = source_recon_results.mni_coord;
          results.mask_indices_in_source_recon = 1:length(results.mni_coord);
        end
    end;
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% setup trial list

D=source_recon_results.BF.data.D;
triallist=[];
for i=1:length(source_recon_results.source_recon.conditions), % indexes conditions/triggers within subset

    trigname=source_recon_results.source_recon.conditions{i};
    Ntrialspercond=length(D.pickconditions(trigname)); %% number of trials for this condition
    triallist=[triallist , D.pickconditions(trigname)];           

end;
Ntrials   = length(triallist);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% setup time indices

if exist('noreject','var') && noreject == 1 
  source_recon_time_indices=1:length(source_recon_results.samples2use);
else
  source_recon_time_indices=find(source_recon_results.samples2use);
end
source_recon_times = D.time(source_recon_time_indices);
results.times=source_recon_times;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% setup channels

modality='MEG';
chanindmeg = strmatch(modality, D.chantype);
chanind = setdiff(chanindmeg, D.badchannels);
if isempty(chanind)
    error(['No good ' modality ' channels were found.']);
end

D=osl_get_oat_sensordata(source_recon_results);

%%%%%%%%%%%%%%%
%% load in sensor data
results.D_sensor_data=D;
results.chanind=chanind;

sensor_data=D(chanind, source_recon_time_indices, triallist);

NK=numel(source_recon_results.BF.inverse.W.MEG); 

if NK>1
    classchanind=find(strcmp(D.chanlabels,'Class'));
    if(isempty(classchanind)),
        error(['No ''CLASS'' chanlabel in: ' D.fname]);
    else
        for kk=1:NK,        
            class_samples_inds{kk} = (D(classchanind,source_recon_time_indices,triallist)==kk);
        end;   
    end;
else       
    class_samples_inds{1}= ones(1,length(source_recon_time_indices),length(triallist));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% recon

S2=[];
S2.source_recon_results=source_recon_results;
S2.chanind=chanind;
S2.class_samples_inds=class_samples_inds;
S2.sensor_data=sensor_data;
S2.num_freqs=1;

Ntrials=size(class_samples_inds{1},3);
for kk=1:NK,  
    for tri=1:Ntrials,        
        S2.timeinds{kk,tri}=find(class_samples_inds{kk}(1,:, tri)); % time indices for class kk
        S2.sensor_data_sub{kk,tri}=sensor_data(:,S2.timeinds{kk,tri},tri);
    end;
end;

for vv=1:length(results.mask_indices_in_source_recon),
    S2.voxind=results.mask_indices_in_source_recon(vv);
    
    [ dat wnorms wnorm wnorms_nai wnorm_nai ] = get_voxel_recon_timecourse( S2 );

    if sum(isnan(dat(:)))==0,
      results.source_timecourses(vv,:)=dat;
      results.wnorms(vv,:)=wnorms;
      results.wnorms_nai(vv,:)=wnorms_nai;
    end
end;
    
